clear all

%  load hedonic data
load hed_data_log


dumyear = dummyvar(year);
 
M = max(hed_m);
b1_mm_l = zeros(M,1);  b2_mm_l = zeros(M,1);
vMin = zeros(M,3); vMax = zeros(M,3);
meanxb = zeros(M,1);

hedrestable = zeros(7,M);
hedsetable = zeros(7,M);
Nm = zeros(M,1);
Rsq = zeros(M,1);
residp = zeros(length(logp),1);

for m = 1:M
    index = find(hed_m==m);
    temptract = tract(index);
    dumtract = dummyvar(temptract);
    temp = mean(dumtract);
    index2  = find(temp>0);
    dumtract = dumtract(:,index2);    
    dumtract = dumtract(:,2:end);
    cons = ones(length(index),1);
    [temp ,~,r,~,stats] = regress(logp(index), [  v_crime(index) hv_crime_sq(index) X(index,:) cons dumtract dumyear(index,2:end) ]); %dumyear(index,:)
    residp(index) = r;
    se = sqrt(diag(inv(length(cons))*(r'*r)*inv([  v_crime(index) hv_crime_sq(index) X(index,:) cons dumtract dumyear(index,2:end) ]'*[  v_crime(index) hv_crime_sq(index) X(index,:) cons dumtract dumyear(index,2:end)])));
    hedrestable(:,m) = temp(1:7);
    hedsetable(:,m) = se(1:7);
    Nm(m) = length(r);
    Rsq(m) = stats(1);
    b1_mm_l(m) =temp(1);
    b2_mm_l(m) =temp(2);
    meanxb(m) = mean([X(index,:) cons dumtract  dumyear(index,2:end)])*temp(3:end);
    
    vMin(m,1) = min(v_crime(index));  vMin(m,2) = prctile(v_crime(index),1); vMin(m,3) = prctile(v_crime(index),5);
    vMax(m,1) = max(v_crime(index));  vMax(m,2) = prctile(v_crime(index),99); vMax(m,3) = prctile(v_crime(index),95);
end
fvMin = vMin;
FvMax = vMax;
% Xvars are age lotsize sqft num_rooms
matrix = zeros(12,5);
matrix(1,:) =  1000*hedrestable(1,:);  matrix(2,:) =  1000*hedsetable(1,:); %vc
matrix(3,:) =  1000000*.5*hedrestable(2,:);  matrix(4,:) =  1000000*.5*hedsetable(2,:); % vcsq
matrix(5,:) =  1000*hedrestable(5,:);  matrix(6,:) =  1000*hedsetable(5,:); % sqft
matrix(7,:) =  1000*hedrestable(4,:);  matrix(8,:) =  1000*hedsetable(4,:); % lot
matrix(9,:) =  hedrestable(3,:);  matrix(10,:) = hedsetable(3,:); % age
matrix(11,:) =  hedrestable(6,:); matrix(12,:) = hedsetable(6,:); % num rooms

matrix = matrix;

rowLabels = {'Violent Crime Rate', '', 'Violent Crime Rate Squared', '', 'House Square Footage', '', 'Lot Square Footage', '', 'House Age', '', 'Number of Rooms', ''};
matrix2latex_a2(matrix, 'hedres.tex', 'rowLabels', rowLabels, 'columnLabels', [], 'alignment', 'c', 'format', '%-6.5f', 'size', []);


clear  index temptract dumtract temp index2 bint r rint stats v_crime hv_crime_sq logp year tract X




%  load transition prob data
load tprob_data


M = max(tprob_m);
pxest = zeros(M,3);sigrhoest=zeros(M,1);
pxse = zeros(M,3);

hedrestable = zeros(3,M);
hedsetable = zeros(3,M);
Nmtprob = zeros(M,1);
Rsqtprob = zeros(M,1);

vc = zeros(M,3,3);
residt = zeros(length(v_crime),1);

for m = 1:M
    m;
    index = find(tprob_m==m);
    cons = ones(length(index),1);
    [temp ,~,r,~,stats] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index)]);
    pxest(m,:) = temp;
    residt(index) = r;
    sigrhoest(m) = sqrt(inv(length(index))*(r'*r));
    pxse(m,:) = sqrt(diag(inv(length(cons))*(r'*r)*inv([ones(length(index),1) Lv_crime(index) tprob_t(index)]'*[ones(length(index),1) Lv_crime(index) tprob_t(index)])));
    vc(m,:,:) = inv(length(cons))*(r'*r)*inv([ones(length(index),1) Lv_crime(index) tprob_t(index)]'*[ones(length(index),1) Lv_crime(index) tprob_t(index)]);
    hedrestable(:,m) = temp;
    hedsetable(:,m) = pxse(m,:);
    Nmtprob(m) = length(r);
    Rsqtprob(m) = stats(1);
end
pxest = pxest';
pxse = pxse';


% Xvars are cons lagged crime t
matrix = zeros(6,M);
matrix(1,:) =  hedrestable(1,:);  matrix(2,:) =  hedsetable(1,:); %cons
matrix(3,:) =  hedrestable(2,:);  matrix(4,:) =  hedsetable(2,:); % lagged crime
matrix(5,:) =  hedrestable(3,:);  matrix(6,:) =  hedsetable(3,:); % t

matrix = matrix;

rowLabels = {'Constant', '', 'Lagged Violent Crime Rate', '', 't', ''};
matrix2latex_a2(matrix, 'tprobres.tex', 'rowLabels', rowLabels, 'columnLabels', [], 'alignment', 'c', 'format', '%-6.4f', 'size', []);




yvals = [1:1:19]'; x= 0; 
% with AR(1), xbar is linear in x. x can be evaluated at any value.
% if xbar is non-linear in x. x needs to be evaluated at every value.
bigtheta0 = zeros(5,length(yvals));
for y=1:length(yvals);
gamma = ones(M,1);
dgammadrho = zeros(M,1); %for delta method
beta = 0.95;
v_crime_bar = x*ones(M,1);
extp=x*ones(M,1);
for m=1:M
B = 1;
for t=2:7 
extp(m) =  pxest(1,m) + pxest(2,m)*extp(m)  + pxest(3,m)*(yvals(y)+t-1);
temp = beta^(t-1);
temp2 = (beta^(t-1))*( pxest(2,m)^(t-1));
temp3 = (beta^(t-1))*( (pxest(2,m)^(t-2))*(t-1));
B = B+temp;
v_crime_bar(m) = v_crime_bar(m) + temp*extp(m);
gamma(m) = gamma(m) + temp2;
dgammadrho(m) = dgammadrho(m) + temp3;
end
end
v_crime_bar = v_crime_bar/B;
gamma = gamma/B;
dgammadrho = dgammadrho/B;
theta0=v_crime_bar-x*gamma;
bigtheta0(:,y)= theta0;
end

theta0 = bigtheta0(:,10); 

aaa=diff(bigtheta0')';
theta0con = aaa(:,1);

%%%
phi_res(1,:) = bigtheta0(:,1)';  
phi_res(2,:) = theta0con';

phi_res;

gammase = hedsetable(2,:)'.*dgammadrho;  %delta method.  se(gamma(rho))  = se(rho)*(d gamma/ d rho)
clear  temp  r  t v_crime_bar y x temp2 index extp m


%  load consump data
load consump_data


M = max(con_m);

for m = 1:M
    index = find( con_m==m & (b1_mm_l(m) +  b2_mm_l(m)*a < 0) );
    vMin(m,1) = min(a(index));  vMin(m,2) = prctile(a(index),1); vMin(m,3) = prctile(a(index),5);
    vMax(m,1) = max(a(index));  vMax(m,2) = prctile(a(index),99); vMax(m,3) = prctile(a(index),95);
end

Na = length(a);


save datafromstata

rentd=zeros(Na,M); 
for m = 1:M
    rentd(:,m) = (b1_mm_l(m) + b2_mm_l(m)*a).*price;
end
temp = dummyvar(con_m);
rentd = sum(rentd.*temp,2);
rentd = .075*rentd;
index = find(rentd<0);
con_m = con_m(index);
rentd=rentd(index);
t = t(index);
save est_res_s2 rentd  M con_m t
